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ABSTRACT 

Scattered light images of transition discs in the near-infrared often show non- 
axisymmetric structures in the form of wide-open spiral arms in addition to their 
characteristic low-opacity inner gap region. We study self-gravitating discs and in¬ 
vestigate the influence of gravitational instability on the shape and contrast of spiral 
arms induced by planet-disc interactions. Two-dimensional non-isothermal hydrody- 
namical simulations including viscous heating and a cooling prescription are combined 
with three-dimensional dust continuum radiative transfer models for direct compari¬ 
son to observations. We find that the resulting contrast between the spirals and the 
surrounding disc in scattered light is by far higher for pressure scale height variations, 
i.e. thermal perturbations, than for pure surface density variations. Self-gravity effects 
suppress any vortex modes and tend to reduce the opening angle of planet-induced 
spirals, making them more tightly wound. If the disc is only marginally gravitationally 
stable with a Toomre parameter around unity, an embedded massive planet (planet- 
to-star mass ratio of 10 -2 ) can trigger gravitational instability in the outer disc. The 
spirals created by this instability and the density waves launched by the planet can 
overlap resulting in large-scale, more open spiral arms in the outer disc. The contrast 
of these spirals is well above the detection limit of current telescopes. 
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1 INTRODUCTION 

Young forming protoplanets may leave observational sig¬ 
natures in their parental disc in the form of distinct gaps, 
vortices, warps, or spiral arms. This is thought to happen 
during the final stages of protoplanetary disc evolution. 
There is a peculiar group called transition discs, whose 
spectral energy distribution (SED) and sub-millimetre 
(sub-mm) observations suggest that the inner disc is 
strongly depleted of dust. Spiral arms have thus far been 
only observed in this kind of discs, suggesting that the gap 
formation mechanism is related to spiral arm formation. 
This needs, however, to be treated carefully since there 
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might be an observational selection effect. Observations in 
the near-infrared (NIR) and (sub-)mm regime of transition 
discs show a variety of asymmetric features. At sub-mm 
wavelengths asymmetries are seen for example in the 
system Oph IRS 48 (van der Marel et al. 2013), SR 21, HD 
135344B/SAO 206462 (Perez et al. 2014), and HD 142527 
(Casassus et al. 2013; Fukagawa et al. 2013), for the latter 
a low-mass stellar companion has been found (Biller et al. 
2012; Close et al. 2014). Since discs are optically thick in 
the optical and NIR, high angular resolution scattered light 
images trace small hot dust particles in the upper disc layer. 
Most discs have double spirals with a large opening angle 
that are almost symmetric, e.g. MWC 758 (Grady et al. 
2013; Benisty et al. 2015), LkCa 15 (Thalmann et al. 2014), 
HD 135344B/SAO 206462 (Muto et al. 2012; Garuh et al. 
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2013), HD142527 (Casassus et al. 2012; Rameau et al. 2012; 
Canovas et al. 2013; Avenhaus et al. 2014a) and HD100546 
(Boccaletti et al. 2013; Avenhaus et al. 2014b). 

The origin of the observed spirals is still debated. For 
instance, with a simple analytical description based on 
the spiral density wave theory, the morphology of spirals 
triggered by a hypothetical planet has been studied by 
Muto et al. (2012), Grady et al. (2013) and Benisty et al. 
(2015). The best spiral fitting model requires a high disc 
aspect ratio h = H/r to account for the large pitch 
angles (angle between the spiral arm and the tangent 
circle) implying a very warm disc. Recently, Juhasz et al. 
(2015) focused on the contrast of spirals compared to their 
surrounding background disc. They modelled the spiral 
waves launched by planets by means of locally isothermal 
hydrodynamic simulations as well as analytic descriptions. 
Based on contrast arguments these authors suggested 
that the spiral arms observed are the results of pressure 
scale height perturbations rather than of pure surface 
density perturbations. Dong et al. (2014) also studied the 
observational signatures in transition discs by combining 
two-dimensional two fluid hydrodynamical calculations 
with three-dimensional radiative transfer simulations, but 
rather focused on the observational signatures of gaps 
opened by one or several planets. They stated that density 
waves and streamers inside the planetary gap, produced 
by planet-disc interactions, can be visible in NIR images. 
Besides processes involving planets, it has been examined 
whether other mechanisms can lead to spiral features, 
such as non-ideal magnetohydrodynamical effects (e.g. 
Flock et al. 2015; Lyra et al. 2015) or gravitational insta¬ 
bility (e.g. Lodato & Rice 2004, 2005; Rice et al. 2004). 
Structures induced by the onset of gravitational instabilities 
tend to produce spirals with higher than m=2 azimuthal 
wave number (e.g. Cossins et al. 2009; Forgan et al. 2011). 
Planets drive, depending on the planet-to-star mass ratio, 
a m=l or sometimes even a m=2 mode, but the secondary 
spiral is substantially weaker. Our idea is to combine both 
spiral formation scenarios in order to see if this strengthens 
their amplitude and to produce various spiral morphologies. 

In general, the precondition for a disc to become un¬ 
stable is that the amount of gravitational potential energy 
overcomes pressure and rotational kinetic support. The ax- 
isymmetric stability of a thin disc is determined by the so- 
called Toomre parameter Q (Toomre 1964). The Toomre 
criterion for a disc to be unstable is defined as 


where c s is the sound speed, E g is the gas surface den¬ 
sity and k corresponds to the epicyclic frequency, which 
is equal to the angular velocity f1 = (GM*/r 3 )'^ in 
the case of a Keplerian disc. Increasing the disc mass 
and/or decreasing the disc temperature leads to a lower 
Q. Hence, the disc mass has to be a significant fraction 
of the stellar mass for the disc to be gravitationally unstable. 

Current observational methods to estimate the disc 
mass still have systematic uncertainties. The total disc 


mass is dominated by the gas, which mostly consists of 
H 2 . Since its emission is hard to detect, there is only 
one object (TW Hya) for which hydrogen was directly 
observed in the form of HD (Bergin et al. 2013). Two 
prominent proxies, dust continuum emission and CO 
and its isotopologues, are usually used to estimate disc 
gas masses. The first method is based on observations 
of the mm continuum emission from dust grains (e.g. 
Beckwith et al. 1990; Andrews & Williams 2005), assuming 
a certain dust opacity and gas-to-dust ratio. The generally 
assumed value of 100 for the gas-to-dust ratio may not be 
accurate for all discs, since it strongly depends on the disc 
evolution (e.g. Brauer et al. 2007; Birnstiel et al. 2010). 
Also, from observations it is known that mm-sized dust 
grains and the gas do not necessarily have the same spatial 
distribution in discs (e.g. de Gregorio-Monsalvo et al. 2013; 
Walsh et al. 2014). Furthermore, in this calculation dust 
opacity values at the observed frequency are assumed, 
which may have quite uncertain values in protoplanetary 
discs due to the various grain sizes, compositions, internal 
structures and the presence of ice mantles (e.g. Pollack et al. 
1994; Henning & Stognienko 1996; Semenov et al. 2003; 
Demyk et al. 2013). The second method infers the gas mass 
independently of the dust content using molecular lines 
observations. Due to its high abundance and line strength 
CO and its isotopologues are the most frequently used 
tracers of gas in protoplanetary discs (e.g. Dutrey et al. 
1996; Williams & Best 2014). Uncertainties concern the 
H 2 -CO abundance and the topologist ratios because of 
photo-dissociation and freeze-out processes. Miotello et al. 
(2014) showed that the disc mass may be underestimated 
by up to two orders of magnitude if only a single CO 
isotopologue line is observed and isotope selective effects 
are not properly taken into account. As mentioned above, 
the only direct measurement of hydrogen in discs from 
TW Hya actually suggests that the discs might be more 
massive than previously thought. Therefore, current disc 
mass estimations from gas and dust observations only give 
a very rough approximation and discs around Class II stars 
might be sufficiently massive to be gravitationally unstable, 
as it is expected for Class 0 and I objects. 

All disc processes including planet-disc interactions and 
gravitational instability strongly depend on the cooling effi¬ 
ciency, and therefore, the disc temperature. The latter again 
has an important effect on the shape and contrast of non- 
axisymmetric disc features. However, many of the numerical 
simulations so far are limited due to the assumption of a 
locally isothermal disc. In this paper, a two-dimensional hy¬ 
drodynamical code which includes an energy equation ac¬ 
counting for disc viscous heating and cooling is used to 
simulate planet-disc interactions. Gravitational instability is 
known to produce multi-armed spirals in discs, while planet- 
disc interactions mostly result in one- or two-armed spirals. 
In this work, we study the role of heating and cooling, and 
gravitational instability on the shape and contrast of spi¬ 
ral arms induced by planet-disc interactions. We want to 
investigate what kind of structures the interplay between 
gravitational instability and planet-disc interactions can re¬ 
produce. This combination has not been investigated so far. 
The question we ask here is whether the presence of a planet 
in a marginally gravitationally stable disc can tip it over the 
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limit. The information about the hydrodynamical disc struc¬ 
ture is a prerequisite for our radiative transfer modelling in 
order to link simulation results to observations. The focus is 
on presenting synthetic scattered light images and on inves¬ 
tigating whether they resemble the observations. This paper 
is divided into five parts. In Sect. 2 the modifications to the 
hydrodynamical code and the basic models are introduced. 
The radiative transfer setup is described in Sect. 3. Subse¬ 
quently, the simulation results are presented in Sect. 4. The 
associated discussion and the conclusions are summarized in 
Sect. 5. 


2 HYDRODYNAMICAL SETUP 

The two-dimensional hydrodynamical grid-based code used 
for simulating the planet-disc interactions is based on the 
fargo-adsg version (Baruteau & Masset 2008a,b). This 
modified version of the original FARGO code (Masset 2000) 
implements an energy equation and the disc self-gravity. The 
full treatment of disc heating and cooling is important, since 
in addition to surface density perturbations, temperature 
perturbations occur as a consequence of shocks along the 
spirals. 

2.1 Surface density description 

The gas disc is characterized by an initial gas surface den¬ 
sity profile E g (r), a temperature profile T(r) and a pressure 
profile P(r). The density profile is taken to be a power law 
combined with an exponential cut-off at long radii, specifi¬ 
cally, 



where r c corresponds to a characteristic scaling radius, 
which is set to = 75 au in accordance with results from 
high angular resolution disc imaging in the sub-mm regime 
(Andrews et al. 2010, 2011). In our simulations this is equal 
to 0.3 times the outer boundary of the disc. E giC describes 
the density normalization factor. The surface density index 
S is taken to be 1 in our simulations. With this surface 
density profile it is ensured that most of the disc mass 
(> 60%) is located between ri n and r c , so that artificial 
reflections from the outer grid boundary can be minimized. 

The initial temperature structure is a power law, leading 
for our case of a non-flaring ( H/r = const.) disc to 

T(r) = r - 1 , (3) 

with mean molecular weight fi, proton mass m p , universal 
gas constant 1Z S , gravitational constant G, and where the 
parameter ho defines the initial disc aspect ratio. 

Such a conical disc defines the boundary between a flar¬ 
ing disc, i.e. H/r is a monotonic increasing function of radius 
r, and a self-shadowed disc, i.e. H/r becomes smaller with 
larger radii (Dullemond & Dominik 2004). Considering that 
we intend to analyse scattered light images, in the case of our 


non-flaring, constant opening angle geometry, any pressure 
scale height perturbation with a sufficient amplitude can 
cast a shadow over the remaining outer disc rather easily. In 
this case of grazing incidence an object can cast the largest 
shadow. Furthermore, the background brightness is as dark 
as possible for grazing infall. Therefore, high-contrast spiral 
structures are most easily made in non-flaring discs. Further¬ 
more, Juhasz et al. (2015) showed that modelling a flaring 
disc affects the brightness of the disc in the outer regions, 
but does not improve the visibility of the spirals. From ob¬ 
servations there is also evidence for very low flaring index in 
MWC 758, one of the sources where a two-armed spiral was 
detected. 

2.2 Heating and cooling terms 

The development of gravitational instability leads to a self¬ 
regulation process due to the competition between heating 
from the instability and cooling. After a local temperature 
increase due to a weak shock induced by the spiral density 
waves, the disc cools back down to a predefined background 
temperature on a certain time-scale. It takes at least a dy¬ 
namical time-scale (fdyn ~ fi _1 ) for the disc to hydrody- 
namically react to any increased temperature. The form of 
the energy equation implemented in FARGO is 

^ +V-(ev) = -PX7-v + Q+-Q (4) 

where e denotes the thermal energy per unit area, v is 
the flow velocity, P is the vertically integrated pressure, 
and Q+/Q_ corresponds to the vertically integrated heat¬ 
ing/cooling term. The disc heating is assumed to be due to 
the disc viscosity. For the cooling source term we assume 

Q- = r~ ( 5 ) 

t'cool 

with a cooling time-scale of (cf. Gammie 2001) 

tcool=P^~ 1 , (6) 

where the /3-factor is taken to be a constant. The powerlaw 
background temperature profile to which the disc cools 
down after each viscous heating event is set by the initial 
locally isothermal temperature (cf. Eq. 3). It is assumed to 
be due to a balance between viscous heating and irradiation 
heating on the one hand, and radiative cooling on the other 
hand. With our choice of the background temperature we 
expect that local heating by shocks becomes relevant. This 
shock heating may increase the scale height of the disc, 
which causes local bumps on the disc surface (radial t = 1 
surface with r being the optical depth). This is different for 
a flaring irradiated disc (T(r) oc r -1 ^ 2 ) around a Herbig 
star, where the heating of the disc is dominated by the star 
apart from the very inner disc. It should be pointed out that 
in a flaring disc model a different background temperature 
profile is probably needed. 

The question is how fast the disc is able to cool down 
after a shock. If the cooling is very efficient, i.e. t CO oi <S fi _1 , 
no shock effects will occur and the disc can be assumed to 
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be locally isothermal. For a rough estimation of f CO oi in our 
models, i.e. /3, we do the following calculation. The cooling 
rate of an accretion disc strongly depends on the optical 
depth r, but since we are in the range r 1 , it is determined 
by 


Q-(r) = 2aT* s (r)=2a-T^ d (r), (7) 

T 

where a is the Stefan-Boltzmann constant and T m id corre¬ 
sponds to the mid-plane temperature. The optical depth can 
be calculated by 

T ~ is d K d , (8) 

whereas a dust-to-gas ratio of 0.01 is assumed and for the 
dust opacity the Planck mean value is taken. Further¬ 
more, the vertically integrated thermal energy is determined 
by 


E = 


_J_ kWT_ 

7 — 1 pm p ’ 


(9) 


with Boltzmann constant fen, proton mass m p , adiabatic 
index 7 = and mean molecular weight p, = 2.3. Using Eq. 
5 and dividing f coo i by the orbital time-scale t OI h = 2 7 rfU 1 
yields the /3-factor. Using typical values for the gas density 
(E g (50au) ~ 50gem -2 ) and temperature (T m id(50au) 
~ 12 K), and setting Kd(T mi d) ~ 10cm 2 g -1 gives /3 ~ 5. 


There is a critical cooling time-scale, t CO oi,c, and a cor¬ 
responding critical value of /3, j3 c , below which the disc 
can fragment and form gravitationally bound clumps, rather 
than reaching a steady, gravitoturbulent state. In numerical 
experiments Gammie (2001) showed that disc fragmentation 
can happen for a typical time-scale of t C ooi,c 3 I2 _1 . Using 
three-dimensional smoothed particle hydrodynamic (SPH) 
simulations Rice et al. (2003) generally confirmed this cool¬ 
ing time for fragmentation, even though the /3 C value is up 
to a factor of two higher for more massive discs. However, 
Paardekooper (2012) showed that disc fragmentation is a 
stochastic process and observed it for cooling times up to 
20 fl -1 . In principle, fragmentation is even possible up to 
/3 = 50, but it becomes very rare for such high values. A de¬ 
tailed study on the convergence of the critical cooling time- 
scale with resolution for SPH and grid-based hydrodynamics 
simulations was done by Meru & Bate (2012). They showed 
that reducing the dissipation from the numerical viscosity 
leads to larger values of /3 C at a given resolution. In the fol¬ 
lowing we decided to use /3-values of 1 and 10, for which 
none of our disc models is found to fragment for the se¬ 
lected resolution. This allows the disc to settle down into a 
quasi-steady, self-gravitating state. We note here that, fol¬ 
lowing the argument from Meru & Bate (2012), repeating 
our simulations with a higher resolution could minimize the 
artificial viscosity and, therefore, possibly fragment the disc 
into bound objects. 


2.3 Parameter choice 

The basic model consists of a viscous, non-flaring disc with 
an embedded giant planet. An overview of all the param¬ 
eters can be found in Table 1. For the complete series of 


models, the disc is tapered as described in Sect. 2.1. FARGO 
uses dimensionless units, therefore, the fixed orbital radius 
of the embedded planet of r p = 1 is used as the length scale. 
The disc is divided into 1280 azimuthal and 720 radial grid 
zones, ranging from 0.2 to 10.0. This corresponds to a phys¬ 
ical disc extension from 5 au to 250 au, assuming that the 
planet’s orbital radius is 25 au. The high resolution for the 
hydrodynamical part serves to avoid strong numerical diffu¬ 
sion effects. The outer boundary of the computational grid 
is therefore far enough away from the regions around the 
gap on which we focus our study. Thus, together with the 
choice of a tapered gas surface density, artificial boundary 
condition effects are minimized. The inner boundary condi¬ 
tions are open, allowing for mass outflow at the inner edge. 
The disc is considered to be non-flaring with a constant as¬ 
pect ratio of h = 0.05. Furthermore, for the disc viscosity 
the a-type viscosity (Shakura & Sunyaev 1973) is used with 
a = [10 - 3 ,10 -2 ]. The planet’s potential itself is softened 
over a length that scales with the disc thickness with a fac¬ 
tor of e = 0.6. The time-scale over which the planet mass is 
switched on for the potential evaluation at the beginning of 
each simulation is set to 100 planetary orbits. The results of 
all hydrodynamical simulations are described in Sect. 4.1. 


3 RADIATIVE TRANSFER SETUP 

The output of the two-dimensional hydrodynamical simu¬ 
lations previously presented in Sect. 2, i.e. the gas surface 
density distribution E g (r,p) and the temperature distribu¬ 
tion T(r,p), are taken for further processing in the context 
of three-dimensional radiative transfer modelling. The radia¬ 
tive transfer code RADMC-3D 1 is used to calculate synthetic 
scattered light images in the NIR. For these calculations the 
same aspect ratio, flaring index, and radial and azimuthal 
grid extensions as in the hydrodynamical simulations are 
adopted. In order to avoid low photon statistics, the out¬ 
puts from the hydrodynamical simulations are interpolated 
to a grid with a lower resolution (N r x N v = 384 x 512 
grid cells) for the radiative transfer calculations. This does 
not affect the shape of any of the structures obtained in the 
FARGO simulations. In addition, Ng = 128 cells for the polar 
direction are taken. Assuming the vertical density profile to 
be Gaussian and adopting the canonical dust-to-gas ratio for 
the interstellar medium of 0 . 01 , the dust volume density is 
given by 


p(R,ip,z) = 0.01 


£g(_B, 1 p) 

V2nH{R,p) 


exp 


2 H\R,v) 


( 10 ) 


where the spherical coordinates R and z can be converted 
into cylindrical formulas via R = r sin($) and z = r cos(#), 
where 9 is the polar angle. The temperature output after a 
specific orbital evolution from the hydrodynamical simula¬ 
tions is used to determine the pressure scale height H(R, p), 




ks 


p m p G 


T(R,p)R 3 oc Vf. (11) 


1 http://www.ita.uni-heidelberg.de/~dullcmond/software/radmc- 
3d/ 
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The radiative transfer calculations start with comput¬ 
ing the dust temperature structure by means of a thermal 
Monte Carlo simulation using 10 7 photon packages. We 
want to calculate an equilibrium dust temperature consid¬ 
ering the star as the source of luminosity. This temperature 
is the result of a balance between radiative absorption and 
re-emission, i.e. a dust grain acquires as much energy as it 
radiates. This implies that the gas temperature from the 
hydrodynamical simulations is assumed to only influence 
the disc scale height (cf. Eq. 11), but is not explicitly related 
to the dust temperature itself. Hence, a self-enhancement 
effect, such that the further irradiation of a disc bump 
actually influences the scale height at those positions, is not 
included in our study. The main inputs for the radiative 
transfer modelling are the dust density structure from Eq. 
10, dust opacities, and the radiation source. The latter 
is assumed to have typical stellar parameters of a Herbig 
Ae star ( T e g = 9500K, M* = 2.0M 0 , R* = 2.5R 0 ). For 
simplicity, a black body radiation field is assumed. The 
dust in our models consists of silicate grains with a fixed 
grain radius of 0.1 /im. The optical constants are obtained 
from the Jena database based on the work by Jaeger et al. 
(1994) and Dorschner et al. (1995). Scattering is considered 
to be anisotropic with a full treatment of polarization. 
The scattering phase function depends on the scattering 
angle and on the polarization state of the input radiation, 
i.e. the parameters of the full Stokes vector [I,Q,U,V]. 
Linearized polarized intensity images can be calculated 
with PI = (Q 2 + U 2 ) 1//2 . The dust opacity tables are pro¬ 
duced with the BHMIE code of Bohren & Huffman (1984) 
for calculating scattering and absorption by spheres. The 
images are calculated at 1.65 fim using 10 8 photon packages. 

The signal-to-noise ratio of the theoretical images ob¬ 
tained from radmc-3d is limited only by photon statistics. 
In order to simulate synthetic observations from these mod¬ 
els that are comparable to observational data, the images 
have to be convolved with the telescope’s point spread func¬ 
tion (PSF). The diffraction limited PSF for a perfect op¬ 
tical system based on circular elements would be an Airy 
pattern, whose central peak can be approximated by a two- 
dimensional Gauss function. As size of the PSF a full width 
at half maximum (FWHM) of 0"04 is chosen, which is rep¬ 
resentative for if-band observations with the SPHERE in¬ 
strument at the Very Large Telescope (VLT). 


4 RESULTS 

4.1 Hydrodynamical results 

The numerical hydrodynamic simulations include a treat¬ 
ment of heating and cooling, providing a more realistic pic¬ 
ture of the disc gas density and temperature distribution in 
self-gravitating discs. Hence, it is possible to study the effect 
of surface density and pressure scale height perturbations on 
the spiral morphology and contrast in scattered light images. 

4-1.1 Energy equation vs. locally isothermal structure 

For our reference model 1 the simulation considers a non 
self-gravitating disc with a mass of 0.15 M*, a planet-to-star 


Table 1 . Overview of the models and parameters we studied. 


# 

M diB c[M*] 

a 

Mp/M* 

p 

r c [r p ] 

SG 

1 

0.15 

10 -3 

10 -3 

10 

3.0 

no 

2 

0.15 

10“ 3 

10“ 3 

10 

3.0 

yes 

3 

0.15 

10 -3 

10 -3 

1 

3.0 

yes 

4 

0.15 

10“ 2 

10“ 3 

10 

3.0 

yes 

5 

0.15 

10 -3 

10 -2 

10 

3.0 

yes 

6 

0.2 

10“ 3 

10 -3 

10 

3.0 

no 

7 

0.2 

10 -3 

10 -3 

10 

3.0 

yes 

8 

0.2 

10“ 3 

10“ 2 

10 

3.0 

yes 

9 

0.2 

10“ 3 

10“ 3 

10 

4.8 

yes 


Notes. When specified the parameters are given in dimen¬ 
sionless FARGO units, i.e. with respect to the planet’s posi¬ 
tion and star mass. SG stands for self-gravity. For all simu¬ 
lations N r X N v = 720 X 1280 grid cells, a radial range from 
r i n — 0.2 r p to r ou t = 10.0 r p , and open boundary conditions 
at the inner edge are used. Model 1 is taken as the reference 
model considering a non self-gravitating disc. 

mass ratio of 10 -3 and an a-viscosity of 10 -3 . The /3-cooling 
factor is set to 10. In Fig. la a snapshot of the gas surface 
density after 1000 planetary orbits is shown. The embedded 
planet clears out a gap in less than 100 orbits. The disc 
viscosity is sufficiently low and the planet-to-star mass ratio 
is high enough, for the angular momentum flux carried by 
the density waves to overcome the inflow due to viscous 
accretion in the disc, leading to the gap opening near 
the planet’s position. The planet is not able to produce 
a full cavity as seen in observations of transitional discs. 
This might be, however, linked to the possibility that dust 
located in the very inner disc close to the star could not be 
detected, for instance due to grain growth. Furthermore, 
distinct spiral density waves are induced by the planet. The 
dominant feature of this model is a vortex at the outer edge 
of the gap that moves with the local Keplerian speed and 
survives until the end of the simulation. The existence of 
this vortex is in accordance with Ataiee et al. (2013), who 
showed that vortices can be long-lived for moderate disc 
viscosity values and massive planets (cf. also Zhu & Stone 
2014; Fu et al. 2014). The vortex influences the surface den¬ 
sity perturbations induced by the planet close to its position. 

Assuming a locally isothermal equation of state corre¬ 
sponds to an infinitely short cooling time-scale. This means 
that the temperature always remains constant after a den¬ 
sity increase since compressive work is immediately radiated 
away. Therefore, the decisive advantage of the energy equa¬ 
tion simulations over locally isothermal models is that the 
information of the temperature distribution is tracked. In 
Fig. le-h it can be seen that the spirals themselves are hot¬ 
ter than the background disc. The temperature of the spiral 
depends on the cooling time-scale. If the cooling time-scale 
is long, the disc material is unable to cool down to its pre¬ 
shock equilibrium temperature before the next heating event 
occurs with the passing of the next shock wave. This leads to 
an overall higher temperature and disc scale height, also out¬ 
side of the spiral itself. Another interesting finding is that in 
all temperature plots, best recognizable in Fig.lg, the spiral 
close to the planet’s position seems to split up in two parts. 
Outwards of the inner spiral and inwards of the outer spi- 
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Figure 1. Surface density (top row) and temperature (bottom row) maps for a planet-to-star mass ratio of 10“ 3 after 1000 planetary 
orbits for different cooling factors (3 and viscosity values a. The disc mass corresponds to 0.15 M*. Note the difference in the x-axis 
scaling. Panels (a) and (e) show the only model for which self-gravity is turned off. The planet’s position is marked with the white and 
black cross, respectively. 
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Figure 2. Surface density map for a planet-to-star mass ratio of 10 —3 after 1000 planetary orbits for a model without (a) and with 
self-gravity ( b ). Panels (c) and (d) show the corresponding temperature distributions. The disc mass corresponds to 0.2M*, /3 is taken 
to be 10. 



Figure 3. Surface density maps (a,b) for a planet-to-star mass ratio of 10 2 after 1000 planetary orbits. The disc mass corresponds to 
0.15 M* (a) and 0.2 M* (6), respectively. The temperature plots are shown in panels ( c,d ). 
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Figure 4. Disc mass evolution for models with different disc and embedded planet masses (left). Azimuthally averaged Toomre parameter 
Q as a function of radius (middle and right panels). For radii smaller than 20 AU Q has very large values, which are not relevant for 
our study and therefore not plotted here. The disc mass in the middle panel corresponds to 0.15 M*, for the right panel to 0.2M+. The 
vertical dashed line represents the position of the planet. The horizontal dashed line at Q=1 shows the critical Toomre value below which 
the disc is unstable to axisymmetric perturbations. A shearing disc becomes unstable to non-axisymmetric perturbations even at slightly 
larger values around Q ~ 1.5 — 2. 


ral wake an additional high temperature arm is seen. Those 
structures are located within the planetary horseshoe or¬ 
bit leading to a rather circular shape. They are caused by 
shock fronts going through the U-turn of the horseshoe or¬ 
bit. We can check whether a planet can generate a strong 
enough perturbation in the temperature to meet the crite¬ 
ria for spiral detection obtained from analytical modelling by 
Juhasz et al. (2015). They found that a relative change of at 
least 0.2 in pressure scale height ( SH/H ) is required for the 
spirals to create detectable observational signatures. Calcu¬ 
lating this relative change for our simulation with /3 = 10 
gives 0.2 at the ip = 270° axis. This value increases to a max¬ 
imum of 0.6 near the planet at a zero azimuth angle. Hence, 
we expect to see the spirals in scattered light (cf. Sect. 4.2). 


4-1-2 Effect of self-gravity 

In Fig. 1 the gas surface density ( a-d ) and temperature 
(e-h) maps of four different models for a disc mass of 
0.15 M* are displayed. The leftmost panel represents the 
reference model with a planet-to-star mass ratio of 10~ 3 
after 1000 planetary orbits and a cooling time of 10 fi -1 . 
The other panels additionally include self-gravity and only 
differ in the /3-factor and a-viscosity value. Models 2 and 
3 (Figs. 16, c) are quite different from the reference model 
without self-gravity. The vortex present in the latter case 
is smeared out very quickly due to self-gravity effects and 
no interaction between the vortex and density waves can 
occur. Additionally, the gap around the planet’s orbit is 
much more distinct. In contrast to very low-mass planets 
driving a single one-armed spiral (Ogilvie & Lubow 2002), 
two arms that spiral outward of the planet’s position with 
a mutual azimuthal shift are present. Note that this m=2 
mode is not a result of self-gravity, but caused by the 
massive planet itself. The primary spiral arm launched 
at the planet’s position is much stronger than the one 
azimuthally shifted which has a quite low density contrast. 
Furthermore, it is noticeable that the self-gravity models 


show more twisted spirals with lower pitch angles. This can 
be explained by the fact that the disc self-gravity modifies 
the positions of the Lindblad resonances, between which 
density waves can propagate. It shifts the location of the 
effective resonances closer to the planet (Pierens & Hure 
2005). 

To study the disc stability we show the azimuthally 
averaged Toomre parameter Q as a function of radial dis¬ 
tance in Fig. 4. The middle panel only includes the Toomre 
profiles for models with a disc mass of 0.15 M*. We com¬ 
pare the models with an embedded planet to the case of a 
self-gravitating disc without a planet (solid black vs. dot¬ 
ted red line). It can be seen that initially Q is always larger 
than 1 for all models, thus the disc by itself would not be¬ 
come gravitationally unstable. Adding a sufficiently massive 
planet considerably reduces Q and puts the disc just at the 
limit of being gravitationally unstable between 50 and 100 au 
from the star (blue dotted line). The explanation for this 
behaviour is the following. When a planet is present in the 
disc it will carve a gap around its orbit decreasing the sur¬ 
face density in the gap. Since most of the material from the 
gap is pushed away due to angular momentum exchange be¬ 
tween the planet and the disc the surface density of the disc 
will increase outside of the gap edge. Hence, a sharp density 
jump is created. Since Q is proportional to SA 1 , it is strongly 
increased near the planet’s position within the gap, but can 
be reduced just behind the gap edge, where overdensities 
grow. Thus, if Q is low enough at the outer edge of the gap 
the local increase of surface density due to gap formation 
might tip the balance and lower Q below unity. The strong 
increase of Q from ~ 100 au on is caused by the exponential 
mass tapering and the disc mass loss (Fig. 4, left panel). The 
issue for a disc mass of 0.15 M* is, however, that Q is not low 
enough at larger radii in order to let self-gravity affect the 
disc dynamics in form of additional gravitational instability 
structures. How a larger disc mass or an increased critical 
taper radius affect these results is discussed in Sect. 4.1.5. 
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4-1-3 Effect of cooling time-scale 




model #9: M rtlaW = 0.2 M* 


According to Eq. 6, for constant /?, t coo i is proportional to 
r 3 ' 2 , so that the cooling is faster in the inner disc. This means 
that the temperature profile returns back to the locally 
isothermal background values faster in these inner parts, but 
in any case not faster than the local dynamical time-scale. 
In general, we expect to see stronger temperature changes 
for longer cooling times, i.e. larger /3-values. The effect of 
the cooling time-scale can be seen by comparing models 2 
(/3 = 10) and 3 (/3 — 1). A very fast cooling (/3 < 1) is very 
close to the locally isothermal case. In Fig. 1 c,d the final sur¬ 
face density and temperature show lower values in absolute 
strength and contrast along the spirals compared to model 
2 ( b,f ). This slight drop in temperature also decreases the 
pitch angle, leading to more tightly wound spirals which are 
smeared out at larger radii. A larger /3 value implies less 
effective cooling, which makes the whole disc, i.e. also the 
spirals, hotter. The higher the temperature, the faster the 
waves spread in the disc, since the sound speed is propor¬ 
tional to the square root of the temperature. Hence, more 
mass flows through the inner disc boundary for /3 = 10 (see 
Fig. 4, left panel). 

4-1-4 Effect of viscosity 

Another important parameter to explore is the viscosity, 
for which the gap profile dependence is already known (e.g. 
Crida et al. 2006). If the viscosity decreases the gap becomes 
gradually deeper. However, the dependence of the gap char¬ 
acteristic on viscosity is less sensitive in numerical simula¬ 
tions compared to theoretical calculations. In simulations 
there is no simple balance between viscous, pressure and 
gravitational torque, since parts of the latter are transported 
away by density waves (Papaloizou & Lin 1984; Rafikov 
2002). Nevertheless, by comparing our model 2 (a — 10 -3 ) 
with 4 (a = 10” 2 ) in Fig. 1, one can see that as viscosity 
increases the gap is filled with more gas. Since the primary 
spiral is launched at the planet’s position within the gap 
this also influences the spiral contrast, which is consequently 
reduced. Furthermore, a higher a-viscosity smears out the 
spiral structure and reduces the absolute surface density at 
the spiral position. Therefore, we use moderate values of 
viscosity for the rest of our models (a = 10“ 3 ). 

4-1-5 Effects of disc and planet mass 

For the simulations 6 and 7 presented in Fig. 2 the disc 
mass is increased to 0.2 M+, but the planet-to-star mass 
ratio of 10 -3 is kept. The same trends as previously 
described for the lower disc mass in Sect. 4.1.2 can be seen. 
Self-gravity (model 7) destroys the vortex and enhances 
the strength of the two-armed m=2 spiral considerably. 
This is consistent with results from Lin & Papaloizou 
(2011), who showed that for sufficiently large disc mass 
and therefore sufficiently strong self-gravity the vortex 
modes are suppressed. Instead, new global spiral modes 
develop. Self-gravity effects become even more dominant 
for higher disc masses. The Toomre parameter profile for 
this model is plotted in the right panel of Fig. 4 (solid and 
dashed black lines). The disc is marginally gravitationally 
unstable from the beginning of the simulation between the 
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Figure 5. Surface density map (a) for a planet-to-star mass ratio 
of 10~ 3 after 1000 planetary orbits for the model with a critical 
taper radius of 4.8r p . The disc mass corresponds to 0.2M+. The 
temperature structure is shown in panel (b). 

planet’s position at 25 au and ~60au. This is, however, 
not sufficient to induce gravitational instability, since the 
outer disc parts remain stable. Even the perturber with a 
planet-to-star mass ratio of 10 -3 is not capable to reduce 
Q to a value < 1 further out from 60 au (dashed green 
line). The general effect of increasing the planet-to-star 
mass ratio is illustrated by comparing Fig. 3 to Figs. 1 and 
2. For both disc masses the higher planet mass forces the 
gap to become eccentric and a distinct two-armed spiral is 
not seen anymore. For a 10~ 2 M* planet mass Q is below 
or close to 1 even for the outer disc region between 100 
and 150 au (red dashed line in Fig. 4). The consequence 
can be recognized in Fig. 36. The massive planet is able 
to initiate gravitational instability, visible as a quite open 
spiral arm in the outer disc. Particular attention has to 
be paid to the temperature structure. The massive planet 
produces significant shock heating within the inner 50 au of 
the disc. These temperature changes cause perturbations in 
the vertical structure of the disc and following Eq. 11 this 
strongly influences the characteristics of our NIR scattered 
light images (cf. Sect. 4.2). 

Apart from changing the planet mass the other possi¬ 
bility to trigger gravitational instability in the outer disc is 
to change the distribution of mass as a function of radius 
by increasing the mass in the outer disc for a given disc 
mass. This is done in our simulation with 0.2 M* by shifting 
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Figure 6. Simulated NIR scattered light images in /7-band polarized intensity (A = 1.65 fin 1 ). All models consider a disc mass of 0.15 M*. 
(a) corresponds to the reference model 1 without self-gravity and shows the image at original resolution as calculated with the radiative 
transfer code RADMC-3D. All other images {b-d) are convolved with a Gaussian beam using a FWHM of Cf/04 (at 140 pc distance), which 
is representative for observations with SPHERE/VLT in the //-band. The central 0 / /l of the image were masked to mimic the effect of 
a coronagraph similar to real observations. 
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Figure 7. Simulated NIR, scattered light images in //-band polarized intensity (A = 1.65 fim). All images are based on a 0.2 M* disc. 
They are convolved with a Gaussian PSF with FWHM of 0/04 (at 140 pc distance). For panel (c) the critical taper radius is set to 4.8r p . 
The central black region (O'/l) mimics the effect of a coronagraph similar to real observations. 


the critical mass taper radius outwards to 4.8 r p (=120 au). 
The corresponding surface density and temperature plots 
are illustrated in Fig. 5a, b. A more open, global spiral 
pattern is seen. These large scale spirals reflect the onset 
of gravitational instability, but close to the planet they 
overlap with the planetary density waves. This instability 
is supported by the Toomre parameter profile in Fig. 4, 
right panel (blue dashed-dotted line). The disc is initially 
gravitationally unstable in most locations. 

When the disc is at the limit of being gravitationally 
unstable only a slight increase in disc mass changes the 
non-axisymmetric structures considerably. This effect can 
be seen in Fig. 3. For a disc mass of 0.2 M* (Fig. 36) a spiral 
arm in the outer disc induced by gravitational instability is 
present. Its pitch angle is significantly higher than for the 
planetary induced spirals from models with a lower disc mass 
(Fig. 3a). 

4.2 Synthetic radiative transfer images 

We aim to investigate the characteristics of the main ob¬ 
servational features in transition discs, i.e. gaps and spiral 


structures, for which full three-dimensional radiative trans¬ 
fer simulations are required. We are interested in the number 
of spiral arms, their pitch angle, and the brightness contrast 
between the spiral features and the surrounding background 
disc to investigate the observability of spirals. We expect 
that the shock heating along the spiral may increase the 
scale height of the disc locally. This causes bumps on the 
disc surface, which are irradiated by the star, and thus pro¬ 
duce a significant brightness contrast. Based on the model 
setup described in Sect. 3 scattered light images in polarized 
intensity at 1.65 fim are calculated. The disc is nearly face- 
on with an inclination angle of 10°. The wavelength at which 
the following images are made is chosen such that the images 
are directly comparable to //-band observations made with 
SPHERE/VLT and HiCIAO/Subaru. For a realistic com¬ 
parison to observations the radiative transfer images are first 
convolved with a 0V04 Gaussian beam, assuming the source 
to be at 140 pc. Then, each pixel intensity value is multi¬ 
plied with r 2 , where r is its distance from the star, in order 
to compensate for the falloff of the stellar irradiation. This 
way the outer disc structures become much better visible. 
A coronagraph is mimicked by masking the inner 0V1 of the 
disc (14au at 140pc distance). The images are normalized 
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Figure 8. Radial profiles of the normalized polarized surface 
brightness belonging to the calculated images in Figs. 6 and 7, 
respectively. Since the focus is on measuring the contrast of the 
outer spirals, the profiles are only plotted beyond an inner circle 
with a radius of 072. 


to the highest disc surface brightness while using a dynamic 
range of 1000 for plotting the images in logarithmic scale. 


4-2.1 Images from models without self-gravity 

The planetary m=l spiral is clearly visible in the scattered 
light image at infinite resolution based on the reference 
model 1 without self-gravity (see Fig. 6a). The brightness 
contrast between the spiral and the contiguous disc ranges 
from 12 to 30. The large vortex present in the hydrodynami- 
cal simulations is seen in the scattered light as well. This has, 
however, to be interpreted with caution, since it is related to 
the two-dimensional character of our hydrodynamical simu¬ 
lations. Vortices can extend through the whole disc from the 
midplane to the atmosphere (cf. Zhu & Stone 2014), but can 
be only accurately simulated with three-dimensional hydro- 
dynamical models. If the vertical elongation of the vortex 
is small, it is expected to appear less extended in scattered 
light, which traces the disc surface layer. After convolving 
this theoretical image with the Gaussian 0704 PSF the vor¬ 
tex as well as the planet-induced spirals are significantly 
smeared out, but still clearly visible. As can be seen in Fig. 
8a, the contrast of the spirals with the surrounding disc is 
reduced by a factor of about three, giving contrast values be¬ 
tween 3 and 10 (magenta line) for a resolution of 0 / /04. We 
note that the contrast is expected to be further reduced in 
observations due to instrument related and additional noise 
effects. 


f.2.2 Images from models including self-gravity 

In all simulations shown in Figs. 6c, d and 7 a-c self- 
gravitating discs are considered. First, we study the effect 
of self-gravity by comparing the results of models 1 and 
2. As already explained in Sect. 4.1.2, the inclusion of 
self-gravity destroys the vortex immediately. This leads to a 
more distinct primary spiral launched close to the planet’s 
position (see Fig. 6c). Furthermore, the radial brightness 
cross-sections in Fig. 8a illustrate that the clear wiggle 
structure from the models without self-gravity (magenta 
line) is smeared out (green line). Self-gravity forces the 
spiral close to the planet to be more tightly wound. This 
in turn reduces the total number of individual spiral 
windings seen in the scattered light, since they cannot be 
distinguished anymore. The contrast of the spiral with the 
surrounding disc varies between 1.5 to 10 for a resolution of 
0704. 

The effect of increasing the planet-to-star mass ratio 
is displayed in Figs. 6 d and 7b. The absolute brightness 
near the planet is increased, leading to several bright spots. 
Moreover, the gap is not fully detectable anymore in the 
scattered light. The irregular structure close to the planet’s 
position is caused by the perturbation in temperature (cf. 
Fig. 3c, d). Thus, this heating strongly increases the disc 
scale height within the gap. This elevated region is now 
irradiated by the star hiding the gap in the scattered light 
image. At the same time structures characteristic of grav¬ 
itational instability are visible. They occur as quite open 
spirals reaching also outer disc parts. The spiral contrast 
in the outer disc is extremely high with values of up to 
40. A second possibility to create gravitational instability 
structures is a model with a disc that is gravitationally 
unstable from the beginning, i.e. no massive planet has 
to work as a trigger. The resulting scattered light image 
in Fig. 7c shows a quite symmetric spiral for the whole 
disc. The planet-induced spiral is only visible close to the 
planetary orbit and the outer disc is dominated by the 
spirals triggered by gravitational instability. For this model 
the polarized surface brightness in Fig. 8 b has a clear 
wiggle structure (green line) with a brightness contrast of 
approximately 3. This contrast is similar to the images 
from models without self-gravity, but a factor of 10 smaller 
than for a self-gravitating disc model with a massive planet 
(10 -2 M*). This shows that the brightness contrast between 
the spirals and the contiguous disc is strongly influenced 
by the planet-to-star mass ratio in combination with the 
disc mass distribution. As soon as gravitational instability 
spirals triggered by a planet develop, the contrast increases 
tremendously, the spiral shape is strongly influenced and 
its pitch angle slightly increases. 

Another feature to mention, which is present in all im¬ 
ages, is the shadow of the planet along a position angle of 
PA = —90 deg. This is just an artefact of our simulation 
methods, since we do two-dimensional hydrodynamical, but 
three-dimensional radiative transfer simulations. This means 
that a surface density enhancement, e.g. a circumplanetary 
ring, will increase the total density at all disc heights in the 
radiative transfer models. It is unlikely that even a circum¬ 
planetary disc around a massive planet with several tens 
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of Mj up can cast a shadow over the whole outer disc in its 
full vertical extent. The shadow is expected to be seen in the 
midplane, but the disc atmosphere should remain nearly un¬ 
affected. 


5 DISCUSSION AND SUMMARY 

We present hydrodynamical simulations of planet-disc 
interactions considering viscous heating, cooling, and 
additionally self-gravity effects. Subsequently, the density 
and temperature structure resulting from these simulations 
are used for the follow-up radiative transfer modelling, in 
order to predict synthetic scattered light images from a pro¬ 
toplanetary disc around a typical Herbig Ae star. The focus 
is set on analysing the morphology and detectability of the 
main observational features of transition discs in scattered 
light observations, distinct gaps, and non-axisymmetric 
spiral arms. 

The conclusions of this paper are as follows: 

1. A clear m=2 spiral structure, and in some cases even 
higher azimuthal wave number modes, can be excited 
by a single planet in a self-gravitating disc. The shape 
of the planet-induced spirals depends on the planet 
mass and disc parameters, such as the viscosity, and 
heating/cooling time-scales. 

2. When self-gravity of the disc is taken into account the 
spirals are found to have a smaller pitch angle, making 
them more tightly wound. For sufficiently large disc 
mass, and therefore sufficiently strong self-gravity, the 
vortex modes are suppressed, and the strength of the 
spiral is enhanced considerably. 

3. All our calculated polarized intensity images show 
a planetary gap and spiral arms with a variety of 
morphologies depending on planet mass, disc mass 
distribution, and the influence of the disc’s self-gravity. 
Convolving the images with a circular Gaussian PSF 
with a FWHM of 0'.'04 (typical for current 8-10m class 
telescopes) lowers the contrast by a factor of ~ 3, but 
the spiral features are still visible. 

4. The scale height along the spiral increases due to the 
temperature changes, forming local bumps on the disc 
surface, and thus creates a signal in the reflected light. 
Therefore, the contrast in the images also depends on the 
cooling time-scale via the /3-factor. In a disc where the 
cooling is very efficient (i.e /3 < 1), which approximates 
the case of a locally isothermal disc, any temperature 
increase is immediately suppressed. The remaining 
density perturbation alone is not capable of producing a 
sufficiently strong spiral contrast in scattered light to be 
observable with current state-of-the-art telescopes (cf. 
Juhasz et al. 2015). 

5. Our result of a non-self-gravitating disc with an embed¬ 
ded planet shows that the combination of planet-induced 
density and scale height perturbation along the spiral 
due to shock heating is sufficient to achieve the contrast 
between the spiral and the contiguous disc necessary 


to be detectable in observations. The contrast values 
estimated for a resolution of 0 / . , 04 are in a range from 3 
to 10, which is consistent with the brightness contrast of 
spirals seen in recent observations. 

6. There are two models in our series for which the Toomre 
parameter profiles allow the creation of spirals due to 
gravitational instability. However, their initiation mech¬ 
anism is different in both cases. For model 9 the disc 
itself is massive enough, and the mass distribution allows 
it to be gravitationally unstable from the beginning of 
the simulation. Spirals are visible in the whole disc, a 
difference is, however, seen in the pitch angle between 
the spirals close to the planet and in the outer disc. The 
spiral is more tightly wound near the planet’s position as 
it is driven by the planet, while the spirals in the outer 
disc are created by gravitational instability showing a 
larger pitch angle. A long-scale spiral with quite regular 
windings is visible in the corresponding scattered light 
image with a contrast of about 3 for a resolution of Cf/CM. 

7. The second possibility is a disc that is initially at the edge 
of being gravitationally unstable. As soon as a sufficiently 
massive planet is included, gravitational instability can 
be triggered by this planet (model 8). A spiral arm with a 
huge brightness contrast of ~ 40 and with a higher pitch 
angle than the spirals driven by a planet in a gravita¬ 
tionally stable disc is generated. The strong temperature 
perturbations near the high-mass planet partly hide the 
low-density gap and dominate the irregular structures in 
this disc region. 

To summarize, perturbations in the vertical structure 
of the disc, i.e. scale height perturbations caused by temper¬ 
ature variations due to planet-induced accretion heating or 
local heating by gravitational instability, are able to create 
the necessary spiral-background contrast seen in scattered 
light observations. The presence of a sufficiently massive 
planet embedded in a marginally gravitationally stable disc 
can lead to a variety of spiral morphologies depending on the 
planet and disc mass. It is possible to create a tightly-wound 
spiral close to the planet and a more open spiral arm in the 
outer disc. However, the explanation of the origin of symmet¬ 
ric, open double-armed spirals seen e.g. in SAO 206462 re¬ 
mains challenging. Recently, three-dimensional global hydro- 
dynamical simulations (isothermal and adiabatic, assuming 
a non self-gravitating disc) combined with radiative transfer 
calculations after a short-time of evolution (10-20 local or¬ 
bits) have shown that inner spirals caused by massive plan¬ 
ets (6 Mj U p at ~100au) can be visible at scattered light 
with similar pitch angle, extension and symmetry as obser¬ 
vational results (Dong et al. 2015). However, at longer times 
of evolution when the disk reaches a quasi-steady state it is 
expected that such a planet opens a visible gap in the outer 
disk (100-200 au) contrary to observational results at NIR 
and (sub-)mm. With our model setup there is either one 
dominant primary spiral in the outer disc when the planet 
is working as a trigger for gravitational instability, or higher 
mode spirals with similar contrasts are seen when the disc 
is already initially unstable. A second planet would be re¬ 
quired in order to create two primary spirals of a m=l or 
m=2 structure with nearly the same contrast. However, for 


MNRAS 000, 1-13 (2015) 


12 A. Pohl et al. 


this scenario the two planets would need to have just the 
right mass ratios and be located exactly at the right radial 
locations and azimuthal angles in the disc to cause symmet¬ 
ric spirals. In addition, other mechanisms apart from self¬ 
gravity may play a role for non-axisymmetric structure for¬ 
mation as well. We note that a caveat of our approach is that 
we do not perform fully consistent radiation hydrodynamics 
simulations, e.g. including the consistent release of accre¬ 
tion energy from the planet via radiation, since this is quite 
computational expensive for a parameter study. Long-term 
simulations in three dimensions that include these effects are 
needed to investigate the differences with our current results 
and for the direct comparison with observations. This will 
be addressed in future work. Additionally, we are aware that 
the cooling law from Gammie (2001), used also in most of 
the previous work about self-gravitating discs, is a simplified 
description. It has to be elaborated in order to calculate the 
disc temperature structure more precisely, which influences 
the local scale height changes, and thus also the brightness 
contrast in scattered light. 
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